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Abstract. We present background concepts of the nuclear density functional theory (DFT) 
and applications of the time-dependent DFT with the Skyrme energy functional for nuclear 
response functions. Practical methods for numerical applications of the time-dependent Hartree- 
Fock-Bogoliubov theory (TDHFB) are proposed; finite amplitude method and canonical-basis 
TDHFB. These approaches are briefly reviewed and some numerical applications are shown to 
demonstrate their feasibility. 



1. Introduction 

The nucleus is a quantum object. The nuclear interaction is not strong enough to localize 
the nucleonic wave function, thus, the nuclear matter stays in the liquid phase even at zero 
temperature pQ. This strong quantum nature in nuclei leads to a rich variety of unique 
phenomena. Extensive studies have been made for constructing theoretical models to elucidate 
basic nuclear dynamics behind a variety of nuclear phenomena. Simultaneously, significant 
efforts have been made in the microscopic foundation of those models. 

Although there have been significant developments in the "first-principles" large-scale 
computation, starting from the bare nucleon-nucleon (two-body & three-body) forces, they are 
still limited to light nuclei with small mass numbers, typically A < 10 [2]. In contrast, the 
density functional theory (DFT) is currently a leading theory for describing nuclear properties 
of heavy nuclei [Sid]. It is capable of describing almost all nuclei, including nuclear matter, with 
a single universal energy density functional (EDF). In addition, its strict theoretical foundation 
is given by the basic theorem of the DFT [6] . Since the nucleus is a self-bound system without 
an external potential, the DFT theorem should be modified from its original form. This problem 
was addressed by recent studies El [9]. In this paper, we present basic concepts of the EDF 
in nuclei and feasible methodologies of the time-dependent DFT. 

2. Backgrounds of nuclear energy density functionals 

In this section, basic properties of nuclei and historical developments in nuclear structure 
theory, which leads to the nuclear energy density functional, are briefly reviewed. To simplify 



the discussion here, we consider an infinite uniform nuclear matter neglecting the Coulomb 
interaction. 



2.1. Basic property of nuclear systems: Saturation and independent-particle motion 

2.1.1. Saturation The volume and total binding energy of observed nuclei in nature are 
approximately proportional to the mass number A. In other words, they have an approximately 
constant density po ~ 0.17 fm -3 and a constant binding energy per particle B/A ~ 8 MeV. Thus, 
extrapolating this property to the infinite nuclear matter neglecting the Coulomb interaction, 
the nuclear matter should have an equilibrium state with po ~ 0.17 fm~ 3 and B/A ~ 16 MeV, 
at zero pressure and zero temperature. This property are called "saturation property", that is 
analogous to the liquid. The most famous and successful model based on this liquid picture of 
nuclei is the empirical mass formula of Bethe and Weizsacker [10t [TT] . This formula contains 
the surface and Coulomb terms in addition to the leading term proportional to A, which well 
accounts of the bulk part of the nuclear binding. 

2.1.2. Independent-particle motion in nuclei There are many evidences for the fact that the 
mean-free path of nucleons is larger than the size of nucleus. The great success of the nuclear 
shell model [12] gives such an example, in which nucleons are assumed to move independently 
inside an average one-body potential. The scattering experiments with incident neutrons and 
protons provide more quantified information on the mean-free path. In fact, the mean free 
path depends on the nucleon's energy, and becomes larger for lower energy pQ. Therefore, it is 
natural to assume that the nucleus can be primarily approximated by the independent-particle 
model with an average one-body potential. For the nuclear matter, this approximation leads to 
the degenerate Fermi gas of the same number of protons and neutrons (Z = N = A/2). The 
observed saturation density of po ~ 0.17 fm -3 gives the Fermi momentum, kp ~ 1.36 fm -1 , 
which corresponds to the Fermi energy (the maximum kinetic energy), Tp = kp/2M ~ 40 MeV. 

2.2. Problems of a mean-field picture 

Evidences of the independent-particle motion encourage us to adopt the mean-field picture of 
nuclei. However, it turned out that the mean-field models cannot describe the nuclear saturation 
property. Let us explain this for the uniform nuclear matter with a constant attractive potential 
V < 0. 

The constancy of B/A means that it is equal to the separation energy of nucleons, S. In the 
Fermi-gas model, it is estimated as 



Since the binding energy is B/A ~ 16 MeV, the potential V is about —55 MeV. It should be 
noted that the relatively small separation energy is the consequence of the significant cancellation 
between kinetic and potential energies. In the mean-field theory, the total (binding) energy is 
given by 



where we assume that the average potential results from a two-body interaction. The two kinds 
of expressions for B/A, Eqs. ([I]) and ([2]), lead to Tp ~ — 5V74 ~ 70 MeV, which is different 
from the previously estimated value (~ 40 MeV). Moreover, the negative separation energy 
(Tp + V > 0) contradicts the fact that the nucleus is bound! 

To reconcile the independent-particle motion with the saturation property of the nucleus, 
the nuclear average potential must be state dependent. Allowing the potential Vi depend on 



S « B/A 



(Tf + V). 
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the state i, the potential V should be replaced by that for the highest occupied orbital Vf in 
Eq. (H|), and by its average value (V) in the right-hand side of Eq. ([2]). Then, we obtain the 
following relation: 

V F « (V) + T F /5 + B/A. (3) 

Therefore, the potential Vf is shallower than its average value. 

Weisskopf suggested the momentum-dependent potential V, which can be expressed in terms 
of an effective mass m* [32] : 

k 2 

Vi = 170 + 1/1-+. (4) 

Actually, if the mean-field potential is non-local, it can be expressed by the momentum 
dependence. Equation leads to the effective mass, m* /m = (1 + Ux/TfY 1 . Using Eqs. 
CO), ©, and (H|), we obtain the effective mass given by 

m* (3 5B 1 \ ~ l n A 

+ ^t^> -0.4. (5) 
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Quantitatively, this value disagrees with the experimental data. The empirical values of the 
effective mass vary according to the energy of nucleons, 0.7 < m*/m < 1, however, they are 
almost twice larger than the value in Eq. ([5]). As far as we use a normal two-body interaction, 
this discrepancy should be present in the mean-field calculation with any interaction, because 
Eq. ([5]) is valid in general for a saturated self-bound system. Therefore, the naive mean-field 
models have a fundamental difficulty to describe the nuclear saturation. 



2.3. Nucleon-nucleon interaction (nuclear force) 

To understand the origin of the problem, properties of the nuclear force provides an important 
key. The saturation property of nuclear density reflects a balance between attractive and 
repulsive contributions to nuclear binding energy. One source of such repulsive effects is the 
nucleonic kinetic energy of the Fermi gas. However, its contribution per particle is proportional 
to p 2 / 3 , which is not strong enough to resist against the collapse caused by the attractive force 
between nucleons. Therefore, the nucleonic interaction must contain a repulsive element. Indeed, 
the phase-shift analysis on the nucleon-nucleon scattering at high energy (-Ei a b > 250 MeV) 
reveals a short-range strong repulsive core in the nucleonic force. The radius of the repulsive 
core is approximately c ~ 0.5 fm. This strong repulsive core prevents the nucleons approaching 
closer than the distance c, which produces a strong two-body correlation, p^{r 1,7%) ~ for 
| r*i — r*2 1 < c. The attractive part of the interaction has a longer range, which can be characterized 
by the pion's Compton wave length A,,-, and is significantly weaker than the repulsion. Thus, 
a naive application of the mean-field calculation fails to bind the nucleus, since the mean-field 
approximation cannot take account of such strong two-body correlations. 

At first sight, this seems inconsistent with the experimental observations. As we mentioned 
in Sec 12. 1\ there are many experimental evidences for the independent-particle motion in nuclei. 
We may intuitively understand that it is due to the fact that the nucleonic density is significantly 
smaller than 1/c 3 . Therefore, the collisions by the repulsive core rarely occur and the system 
can be approximately described in terms of the independent-particle motion. Furthermore, the 
effects of the Pauli principle hinder the collisions, since the nucleons cannot be scattered into 
occupied states. Although the repulsive-core collisions are experienced by only a small fraction of 
nucleons (~ poc 3 ), each collision carries a large amount of energy. Therefore, the repulsive core 
provides an important contribution to the total energy and are responsible for the saturation. 

Another important factor for the independent particle motion is the strong quantum nature 
due to the weakness of the attractive part of the nuclear force. The importance of the quantum 



nature can be measured by the magnitude of the zero-point kinetic energy compared to that of 
the interaction. If the attractive part of the nuclear force were much stronger than the unit of 
h 2 /Mc 2 , the quantum effect would disappear and each nucleon would stay at the bottom of the 
interaction potential. Then, the nucleus would crystallize at low temperature. In reality, the 
attraction of the nuclear force is so weak that it barely produces many-nucleon bound states at 
the relatively low density. 



2-4- From Brueckner theory to EDF 

The nuclear matter theory pioneered by Brueckner gives a hint for a solution for the previous 
difficulty to understand the nuclear saturation. The Brueckner theory may provide a first step 
toward the quantitative treatment to understand the saturation property and the independent- 
particle motion in nuclei. Details of the theory can be found at Refs. |144 115j. 

The basic ingredient of the Brueckner theory is a two-body scattering matrix of particle 1 
and 2 inside nucleus caused by the nuclear force v, 

G(u) = * + v J? GH (6) 

u - Q{Ti + T 2 )Q 

where Tj is the kinetic energy of particle i, Q is the Pauli-exclusion operator to restrict the 
intermediate states, and ui is called a starting energy that depends on energies of particle 1 and 
2. This is called G-matrix [16j . The G-matrix renormalizes high- momentum components in 
the bare nuclear force and becomes an effective interaction in nuclei under the independent-pair 
approximation. The G-matrix reflects an underlying structure of the independent many-nucleon 
system through the operator Q and the starting energy to. Inevitably, the G-matrix becomes 
state (structure) dependent. 

Since the short-range singularity is renormalized in the G-matrix, we can calculate the total 
energy in the independent-particle (mean-field) model, analogous to Eq. ([2]). 
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where = e, + €j, defines the self-consistency condition for the Brueckner's single-particle 
energies, and G{j^j = — Gijji- This is called Brueckner-Hartree-Fock (BHF) theory. The 

validity of the BHF theory is measured by the wound integral k = (tp — (f)\ip — (ft), where \cp) is 
an unperturbed two-particle wave function and |V>) is a correlated two-particle wave function 
in nucleus, k is known to be about 15 %. The BHF calculation was successful to describe the 
nuclear saturation, however, could not simultaneously reproduce empirical values of B/A and 
po, known as a problem of the Coester band [T7]- Its applications to finite nuclei also have 
similar problems to reproduce the energy, radius, and density in the ground state. 

A possible solution to these problems was given by Negele [IB]. Starting from a realistic 
G-matrix, using the expressions for the Pauli operator 

{rir 2 \QW 2 ) = {5(n - rj) - p{n ~ r\)} {6(r 2 - r* 2 ) - p{r 2 - 7 2 )} , (8) 

and the average single-particle energy e[p(r)], the local density approximation is introduced 
to expand the off-diagonal density matrix p(r + s/2,f — s/2) with respect to the relative 
coordinate \s\. Then, a short-range part of the G-matrix, which is not fully understood, is 
phenomenologically added to the energy expression to quantitatively fit the saturation property, 
and finally, the total energy is treated variationally. This procedure is called the density matrix 
expansion (DME) [19]. The state dependence of the G-matrix is now replaced by the density 



dependence. The final result for the total energy, for a uniform nuclear matter, is written as a 
function of the neutron and proton densities, p n and p p , and the kinetic densities, r n and t p ; 
E = E[p 

nj Ppi 7~ni This expression can be generalized for finite nuclei as 

E = j dfH(r), H(r) = H[p n (f), p p (f),r n (r),r p (r)], (9) 

H{r) here is a complete analogue of the Hamiltonian density in the Skyrme EDF (without the 
spin-orbit and Coulomb terms) [20]. Thus, the DME with a microscopic G-matrix leads to a 
variational treatment of a simple EDF of the Skyrme type. 

The essential aspect of DME is in its non-trivial density dependence and the variational 
treatment. Now, the expression for the total energy, Eq. ([2]), should be modified due to these 
non-trivial density-dependent terms. This resolves the previous issue, and provides a consistent 
independent-particle description for the nuclear saturation. In nuclear physics, this was often 
interpreted as the density-dependent effective interaction. In this terminology, the variation 
of the total energy with respect to the density contains re-arrangement potential, dV e s[p]/dp, 
which comes from the density dependence of the effective force V^ff [p]. These terms turn out to 
be crucial to obtain the saturation condition. 

In summary, the failure in the mean-field description of nuclei using phenomenological 
effective interactions can be traced back to the missing state (structure) dependence. The 
EDF approaches take into account the state dependence in terms of the non-trivial density 
dependence. 



3. Time-dependent Hartree-Fock-Bogoliubov theory 

For a quantitative description of heavy nuclei in open-shell configurations, it is necessary to 
include the pairing correlations. This can be done by a straightforward extension of the Skyrme 
energy functional simply added by the pairing energy functional. 

E[p, k] = E Skyrme [p] + E pail [p, k] , (10) 

where k is the pair density. The variation of the energy with the constraint on the average 
particle number leads to the Hartree-Fock-Bogoliubov (HFB) equation |15|. I21j. However, here, 
let us start from the time-dependent Hartree-Fock-Bogoliubov (TDHFB) equation, 

i^n(t) = [n(t),n(t)], (ii) 



where T-L(t) is the HFB Hamiltonian usually written in a form [15] 
The generalized density matrix R(t) is given by 

which is Hermitian and idempotent: 1Z 2 = 1Z. Thus, the eigenvalues of 1Z(t) are and 1. The 

eigenvectors = ^y^j^ , which correspond to the eigenvalue 0, are called "unoccupied" 

quasiparticle (qp) orbitals. For every unoccupied qp orbital, there exits a conjugate "occupied" 
(V*(t)\ 

orbital *ff v (t) = I J > whose eigenvalue is 1. 



We define the following matrix of size 2M x M where M is the dimension of the adopted 
single-particle Hilbert space, 

^ Mt) -\(a-M\V u (t)) a = M + l,-..,2M [ ' 

which collectively represents time-dependent unoccupied qp orbitals (y = 1, • • • , M). The matrix 
for the occupied orbitals ^f(t) are defined in the same manner, with ^ u (t) replaced by ^(i). 
Using these matrices, TZ(t) can be written in a simple form as K(t) = *(i)* t (t) = l-*(i)#t(i). 
The orthonormal property of the qp orbitals is given by tft(t)^(t) = ^t(t)^r(t) = i. Note that 
7Z(t) can be regarded as the projection operator onto the occupied space. Substituting these 
expressions into Eq. (jlip . we obtain the TDHFB equation for the qp orbitals as 

K{t)U^--H(t)\y(t) = 0, Q(t)\&-H(t)\*(t) = 0, (15) 



, &t W J w w i dt 

where Q(t) = 1 — 7Z(t) is the projection operator onto the unoccupied space. Therefore, the 
occupied-occupied and unoccupied- unoccupied matrix elements of id/dt — %{t) are arbitrary at 
every instant of time. This is related to the U (M) invariance; TZ(t) is invariant with respect to 
the transformation among the occupied (unoccupied) qp orbitals, ^f' = ^fU = ^U), where 
U is an M x M unitary matrix. 

Now, let us derive the static Hartree-Fock-Bogoliubov (HFB) equation as a quasi-stationary 
solution of Eq. (|lip . This does not correspond to the solution for dlZ/dt = 0, because the 
solution must have a given value of the average particle number; trp = N. When the system 
is in the superfluid phase (k ^ 0), TZ(t) and H(t) do not commute with the matrix associated 

with the particle number, M = (\ ' l^us, we ma y construct the generalized density R{t) 

defined in the frame of reference rotating in the gauge space. 

R{t) = eM^mm exp(-z^t) = (_^ ( ^-a^ fl^) • (16) 
For this R(t), TDHFB equation (HI]) becomes 

i^R(t) = [H(t)-ntf,R(t)], (17) 

Hit) = exp(ipMt)n(t) eM-^t) = (_ A *(|U^ ^f^) ^ 

Namely, the transformation does not change p and h(t), but modifies n(t) and A(i) by 
multiplying the time-dependent phase e 2 *^ 4 . As far as A(t) linearly depends on n(t) and hit) is 
a functional of a product of n(t) and n*(t), H[R(t)] has the functional form identical to ~H\TZ(t)]. 
The stationary state can be defined by the time-independent density in the rotating frame, 
dR/dt = 0, leading to 

[H-pM,R]=0, (19) 

where \x correspond to the chemical potential that should be determined by the condition on 
the particle number. With a proper choice of the U{M) degrees of freedom, the quasiparticle 
orbitals <]?„ become time- independent eigenstates as [H — pN)Q u = E u & u . The collective 
representations, 2M x M matrices and are constructed in the same manner as Eq. (|14j) . 

The HFB state is a quasi-stationary state in which the pair density and potential are 
rotating in the complex plain with a constant angular velocity 2p . This is a collective motion 
corresponding to the Nambu-Goldstone mode associated with the spontaneous breaking of the 
gauge symmetry, called "pair rotation". 



where 



4. Finite amplitude method 

Although there have been significant efforts to develop a numerical solver of the TDHFB 
equation in last decade \22\ I23j. the computation of the full three-dimensional (3D) dynamics 
is still a difficult and challenging task. A possible simplification is given by the small- 
amplitude approximation. It is well-known that this will lead to the quasiparticle random-phase 
approximation (QRPA). Recently, we have proposed a feasible approach for the linear response 
calculation |24] and it has been extended for the QRPA [25] , The method is called "finite 
amplitude method" (FAM). 

The QRPA equations can be written in a form 



(E v + E v , -oj)X uu ,(uj) 
(E v + E v , +u)Y vv ,(u) 



+ {H l 



where 
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Here, Vi(uj) is an external field and H\(uj) is an induced residual field. The most tedious part is 
the calculation of Hi(u) expanded in terms of X(oj) and Y{oS). The FAM calculates the induced 
residual fields as follows 



H[R„] - H[Rq] 
V 



^ H[R V ] - H[Rq] 



where rj is a small real parameter. Here, i?^(w) 



*;(w)*J(w), with 



($ + r ? $y*(w)) 1 



Equivalently, the FAM formula can be written in terms of the qp orbitals as 



H[%*\]-H[$,&] 
V 



- t H[%,ty]-H[*,&] 
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(23) 
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According to these FAM formulae, the 
residual fields can be calculated by using the 
finite difference with the parameter r\. In this 
way, the calculation requires only the HFB 
Hamiltonian H[R]. 

Computer programs of the FAM have been 
developed for spherical nuclei [25] based on 
the hfbrad [26], and for axially deformed 
nuclei [27] based on the hfbtho [28]. The 
FAM in the 3D grid-space representation was 
also achieved |29t [30] for nuclei in the normal 
phase. Here, we show the isoscalar monopole 
strength function in 124 Sn in Fig. [1] The 
ground state is obtained by the hfbrad using 
the SkM* functional and the volume-type 
pairing. The zero-range nature of the pairing 
functional requires us to truncate the model 

space. The truncation was done according to the cutoff for the qp energy, and results with 
different cutoff energies are shown in Fig. [1] as well. The final results are almost independent 
from the choice of the cutoff. 
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Figure 1. Isoscalar monopole strength function 
in 120 Sn. See text for details. 



5. Canonical-basis TDHFB 

In this section, we present another approximate treatment of the TDHFB. Here, we do not take 
the small amplitude limit, instead, assume the diagonal property of the pair potential. In fact, in 
the stationary limit, this corresponds to the the well-known BCS approximation for the ground 
state [15 . Therefore, it can be regarded as the BCS-like approximation in the time-dependent 
treatment. 

At every instant of time, we may identify the canonical basis in which the density is diagonal, 
Pij{t) = pi(t)5ij. The canonical states are always paired with pu(t) = Pjj(t) = pi(t) and 
Kjj(i) = K i(t)5ij- Here, we introduce an assumption that the pair potential is also diagonal in 
the canonical basis, Ajj = — Aj(t)5^ In the stationary limit, this approximation is identical to 
the usual BCS approximation [15J. Then, we end up the following set of equations 



4^*)) = W*) - •*(*)} !*(*)>• l Wt m = {m ~ ^ )} m > (26) 

ij t Pi{t) = «i(t)A?(t) - <(t)A,(i), (27) 

= {Vi(t) + mi*)) + A i(t) {2f>i(t) - 1} • (28) 

Here, rji(t) and r]i{t) are arbitrary real functions to control the gauge degrees of freedom of 
the canonical states. Equations (|2f)]) . (|27p. and ([25]) are invariant with respect to the gauge 
transformation with arbitrary real functions, 6i(t) and &j(t). 

IM^e^Wfa) and |fc> e *i(*)|fc) (29) 
Kk ^ e -W)+flj(0)«. and A k -> e-^+^Ai (30) 



simultaneously with 



rk(t)->m(t) + ^ and m (t) -+ rn(t) + ^. (31) 



Thus, r]i(t) and rfc(t), control time evolution of the phases of |<^>j(£)), \4>j(t)), Ki(t), and Aj(i). 

The canonical-basis TDHFB equations ([26]) . (f27|) . and ([28]) . with a proper gauge choice 
guarantee the following properties [31]: 

(i) Conservation law 

(a) Conservation of the orthonormal property of the canonical states 

(b) Conservation of the average particle number 

(c) Conservation of the average total energy 

(ii) The stationary solution corresponds to the HF+BCS state. 

(iii) In the small-amplitude limit, the Nambu-Goldstone modes correspond to zero-energy 
normal-mode solutions. 

In this study, although the method is applicable to the large amplitude dynamics, it 
is utilized to study the photoreaction of Xe isotopes in the linear regime. Using the 
instantaneous perturbative external El field, we calculate the time evolution in the 3D grid- 
space representation. Then, the Fourier transform of the time-dependent El moment leads to 
the El response function. The smoothing parameter of T = 1 MeV is used. See Ref. |31] for 
numerical details. The calculated El strength distributions in 132_140 Xe are shown in Fig. [2j 



The ground states of these nuclei are calculated to 
be spherical, except for Xe which has a very small 
deformation of (3 = 0.02. The peak energy of the 
giant dipole resonance (GDR) is roughly constant and 
-Egdr ~ 15 MeV for these isotopes. 136 Xe corresponds 
to the neutron magic number N = 82. The property of 
the GDR does not significantly depend on the neutron 
magicity. In the low-energy region below E = 10 MeV, 
we may notice an onset of a small dipole peak beyond 
N = 82, which appear in 138 Xe and increases in 140 Xe. 
This seems to be mainly due to a drastic decrease in 
the neutron separation energy beyond N = 82. A 
systematic study on the low-energy El strength in this 
mass region is currently under progress. 

It should be emphasized that the computational 
cost of the present real-time approach is significantly 
smaller than the normal QRPA calculation based 
on the diagonalization method. For instance, the 
computational cost of the QRPA calculation in Ref. 
|32j is larger than the present one by several orders of 
magnitude, even though the axial symmetry restriction 
was utilized. This is because the canonical-basis 
method utilize the selected canonical states whose 
number is much smaller than the number of qp states. 
In addition, the real-time method is very efficient for 
the present purpose, because the single time evolution 
produces the nuclear response for the entire energy 
region. 

6. Summary 

We presented basic concepts of the density functional 
theory in nuclear physics, based on the historical 
developments in nuclear many-body theories. The 
energy density functional is established as the 
density-dependent effective interactions, which gives a 
consistent mean-field-type description of the nuclear 
saturation property. The time-dependent density 
functional approach is a powerful method to study 
dynamics of the quantum many-body systems. In 
description of heavy open-shell nuclei, the Kohn-Sham 
orbitals should be extended to the Bogoliubov-type 
quasiparticle orbitals, to include the pair correlations. 
This is essentially identical to the TDHFB theory. 
Although the full treatment of the TDHFB is 
still a challenging task, we presented approximated 
treatments; the finite amplitude method and the 
canonical-basis TDHFB method. Both methods serve 
as an efficient computational approach to dynamical 
properties of heavy superfluid nuclei. 
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Figure 2. El strength distributions in 
Xe isotopes calculated with the SkM* 
functional. 
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